docs/examples/pyaudi/Example 1 - PyAudi - Differential Intelligence.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Differential Intelligence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(original by Dario Izzo - extended by Ekin Ozturk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we show the use of `desolver` for numerical integration following the notebook example here [differential intelligence](https://darioizzo.github.io/audi/notebooks/example10.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing Stuff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyaudi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "PyAudi backend is available.\n",
      "Using numpy backend\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "import os\n",
    "import numpy as np\n",
    "\n",
    "os.environ['DES_BACKEND'] = 'numpy'\n",
    "import desolver as de\n",
    "import desolver.backend as D\n",
    "\n",
    "D.set_float_fmt('gdual_double')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Controller representation and “simulator”"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take as an example the task of learning a robotic controller. In neuro evolution (or Evolutionary Robotics), the controller is typically represented by a neural network, but for the purpose of explaining this new learning concept we will use a polinomial representation for the controller. Later, changing the controller into a NN with weights as parameters will not change the essence of what is done here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Definition of the controller in terms of some weights parameters\n",
    "def u(state, weights):\n",
    "    x,v = state\n",
    "    a,b,c,e,f,g = weights\n",
    "    return a + b*x + c*v + e*x*v + f*x**2 + g*v**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Definition of the equation of motion (our physics simulator propagating the system to its next state)\n",
    "def eom(state, weights):\n",
    "    x,v = state\n",
    "    dx = v\n",
    "    dv = u(state, weights)\n",
    "    return (dx, dv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical Integration - Runge-Kutta 45 Cash-Karp Method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Evolutionary Robotics, Euler propagators are commonly used, but we would like to use a higher order integration scheme that is adaptive in order to minimise computation, and increase the accuracy and precision of the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = D.array([D.gdual_double(0.05*(np.random.uniform()-0.5), _, 7) for _ in \"abcefg\"])\n",
    "x = [D.gdual_double(2*(np.random.uniform()-0.5))]\n",
    "v = [D.gdual_double(2*(np.random.uniform()-0.5))]\n",
    "y0 = D.array(x + v, dtype=D.gdual_double)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rhs(t, state, weights, **kwargs):\n",
    "    return D.array(eom(state, weights))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We integrate the system using the Runge-Kutta-Cash-Karp scheme as the numerical integration system with a dense output computed using a piecewise C1 Hermite interpolating spline. \n",
    "\n",
    "This particular interpolator is used as it satisfies not only the state boundary conditions, but also the gradients and is well suited for approximating the solution continuously up to first order in time.\n",
    "\n",
    "Note that the entire integration is done using gduals and thus the solution and the interpolating spline stored in the OdeSystem instance, `pyaudi_integration`, contains all the information about how the state reached by our robot changes when we change the control weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "edf168b4bf7242bab7f7255292aa3879",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=9000000000.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "pyaudi_integration = de.OdeSystem(rhs, y0=y0, dense_output=False, t=(0, 10.), dt=0.1, rtol=1e-12, atol=1e-12, constants=dict(weights=weights))\n",
    "\n",
    "pyaudi_integration.set_method(\"RK45\")\n",
    "pyaudi_integration.integrate(eta=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de3xcdZ3/8ddnZpJMkubaJE2aSxug9MalhXATEChFiiJVUSgiP/itKz9d0NVVWVxcdfXHb3V3WWQX1qULCC4oi4hauQhYEAQpbVoK9H6lTZq2SZs0SZtMLjPf3x+ZljRNb2RmzmTm/Xw88phzznxzvp+h4Z1vvvM9Z8w5h4iIpD6f1wWIiEhiKPBFRNKEAl9EJE0o8EVE0oQCX0QkTQS8LuBwSkpK3MSJE70uQ0RkVFm6dOku51zpcM8lbeBPnDiR+vp6r8sQERlVzGzL4Z7TlI6ISJpQ4IuIpAkFvohImlDgi4ikCQW+iEiaUOCLiKQJBb6ISJqIyTp8M5sD3AP4gQeccz8c8nwN8AhQGG1zu3Pu2Vj0LdIfjrCvJ0x3X5hQX5hQf5hQX2Rguy9MT//Adk9fhJ7+MOGII+Ig4lz0a2DbOYgMes5nRsBvBHxGwO+LPhoZPh9ZGT6yM/xkZ/oPecwLZpCb6cfMvP5PI3KQEQe+mfmB+4DLgEZgiZktcM6tGtTs28ATzrmfmNk04Flg4kj7ltTRF47Q1tVL274+Wvf1Dnx19dLe1UtHqJ/OUF/0cWC7M9RPR/fAY3df2OvyD+H3GfnBAAXZGeRnZxx4LMnNpGRMFiV5WZRGH0vGDBwLZvi9LltSXCxG+GcDG5xzmwDM7HFgLjA48B2QH90uAJpi0K8kub5whJbOHpo7e2juCB30uGtvz/vBvm8g1A8nK+AjL5hBfjBAXjBAfnYGFQVB8rIyyAsGBkbUWQOj62DATzDDTzDDd+Axa9CxrICfgM/wmWE+8Jnhs/2P72+bQcRBfyRCf9gNfEUi9EccfeEIPf0RunsH/qro7g3T1Tvw10RXbzj6y6mP9u4+Orr7Bx5DfWxr62bX3p7DvtaSMZmML8ymsjD7oMeqomxqS3LJzUraC+NllIjFT1Al0DBovxE4Z0ib7wEvmNmXgVxg9nAnMrObgZsBampqYlCaxINzjrauPpr2dLNtTzc7O0Ls7AjR3DEQ7js7QrR09tDa1cvQD1TzGYwdk0XJmCzG5mZSVZRDcW4mRTmZFOdmUJybRVFuBsW5mRTnZFKQk0FWwJuRr9/A7/MT65wN9YXZva+XXdFffLv29rCzo4ft7d00tnWzbmcnL69tJtQXOej7xuVncULJGGpLczmhJJfaklwmleVRVZSNz6fpIzm6WPwoD/eTNvRzE68DHnbO3WVm5wH/bWanOOcO+ol2zs0H5gPU1dXpsxc9EuoLs6M9dCDQm/YMbDe179/vPiSM/D6jZEwm4/KDVBVlc8aEIsrysijLCzIuf+CxLH8g5AP+9F4rEMzwUxkdwR/O4F+qW1u72LxrH5ta9rF5116ee3c7bV19B9rmZvqZXJ7H1Ip8plTkM7U8j8nleeQFMxLxcmQUiUXgNwLVg/arOHTK5vPAHADn3BtmFgRKgOYY9C/HKdQXpqG1i62tXTS0dtHYFg3ztm627Qmxa2/PId9TmpfF+MJsppTnMWtyGeMHTTuMK8hibG4Wfo0yY8bMBv7Kyc3klMqCQ55v29fLpl372NDcyertnaza3sHv3m7isTe3HmhzQkkuM2oKmVldyMyaIiaX55GR5r9s010sAn8JMMnMaoFtwDzgs0PabAUuBR42s6lAEGiJQd8yjEjE0dzZw9ZoqO8P9v0h39x5cKAHM3wH5ounVuQfCPPxhUEqC7MpLwh6Nq0iwyvKzeTM3EzOnFB04Jhzjqb2EGu2d7CqqYO3G9t5dV0LTy3bBgy8F3JqZQFnTCjinNpizqotJl9/BaQVc0MnWT/IScw+CvyYgSWXDznn7jSz7wP1zrkF0ZU5/wWMYWC65zbn3AtHOmddXZ3T7ZEPb29P/0Gj9MHh3tjWTW//+1MuZjC+IJvq4mxqinOoKc6hOvpVU5zD2NxMLSFMUc45Gtu6Wd6wh7e27mF5QxsrtnXQG47gMzi1soBzTxzLeSeM5ayJxXpjOAWY2VLnXN2wz8Ui8OMh3QO/Pxxhe3toYGTetj/Muw8EfOu+3oPa5wUDTBibQ3XR+4G+P9zHF2aTGdCf8jIg1Bdm2dY2Fm3czRubdrO8YQ99YUfAZ8ysKeSSKWXMmlLG5HF5GgiMQgr8JOScY09X36Aw76KhtftAwG9r66Y/8v6/TcBnVBZlHxTm+8O9pjiHghz9aS4fTFdvP8u27OHPG3fxyroWVjZ1AFBZmM3Fk0uZNaWMD51YQnampvVGAwW+R0J9YRrbumloO3gOfWtrN42tXXT2HLwee2xu5qCplvfDvbooh4qCYNqvbpHE2NkR4uU1zby0ppnXNuyiqzdMMMPHJZPLuPK08VwypZScTE39JCsFfpyEI47mztCBkfnW6Oi8ITpa39EROqh9VsB38Ai9OIfqomxqolMxmj+VZNPTH+bNTa28uGonz63Ywa69PWRn+Jk1tYwrT63g4sllGvknGQX+B9QZ6juwBn3/+vOm6Lr0/RccDZ52MYOK/CBVg+bP979RWl2UQ2leluZEZdQKRxyLN7fyzLtN/H7FDnbt7SUn08+cU8q5pq6ac2qL9fOdBNIr8B97DO64A7ZuhZoauPNOuP76g5o45+gI9dPS+f7VoU3t74f5/oDvHHIJfMBnlBcEB132HqSiIPvASL2yKFvLFyUt9IcjLN7cyu/eaeLpt7fT2dNPTXEOnzmziqvPrGL8ES4qk/hKn8B/7DHczTdjXV0HDvVlBVlwy/dYeMbsA+He3Bk65EpRgMKcDMYX7L+gKDhoPfpAwJfm6eIikaG6e8P8fuV2flnfyJ837sYMLjiphM+eXcNl08bpvacES5/AnzgRtmw55HBTQRmf+/v/oSwvi3H5wQOX/JflZ1Ea3a4oCGoOXWSEGlq7+OXSRp6sb6CpPcT4giA3nDeR686upjAn0+vy0kL6BL7PxyF364KByfXIoSN6EYmPcMSxcPVOHv7ze/x5426CGT4+ObOSmz5Uy+TyPK/LS2lHCvzUGtLW1Aw7wkd33hRJKL/P+Mj0cj4yvZw1Ozp45M/v8dSybfxicQMXTirhlktO4twTxnpdZtpJrcm1O++EnJyDj+XkDBwXEU9MKc/nHz91Gou+dSm3zZnM6u2dzJu/iGv+8w1eXddCss4ypKLUmtKBY1qlIyLeCfWF+Z8lDfznKxvZ3h7i9KoCbp01idlTy7SsMwbSZw5fREaN3v4ITy1r5D/+uJGtrV1MH5/P7VdM4cJJpV6XNqodKfBTa0pHREaNzICPeWfX8NLXL+Kuz5zOnq4+bnhwMTc8+CYrtrV7XV5KUuCLiKcCfh9Xn1nFS9+4iG9/bCrvbmvnyn9/ja8+/hYNrV1HP4EcM03piEhSae/u4/5XNvLga5txDj5/YS1fnnWSbth2jDSlIyKjRkF2BrfNmcIr37yEj58+np/8cSOz73qF36/YrhU9I6TAF5GkVF4Q5K5rTueXXzyP/OwMvvjoMm766RI279rndWmjlgJfRJLaWROLefrLF/CdK6exdEsbl9/9Kne/uO6gj/GUY6PAF5GkF/D7+IsLannp6xdxxanl3LNwPVfd+5pW8xwnBb6IjBpl+UHumTeTB2+sY/e+Xj5x3+sa7R8HBb6IjDqXTh3Hi1/7MB8/fTz3LFzPJ+57nVXRz+KVw1Pgi8ioVJiTyd3XzuD+G86kuTPE3Pte44E/bdJKniOISeCb2RwzW2tmG8zs9sO0ucbMVpnZSjP7eSz6FRG5fHo5L37tIi6ZXMb/fWY1X/hZPXu6er0uKymNOPDNzA/cB1wBTAOuM7NpQ9pMAr4FnO+cmw58daT9iojsV5Sbyf03nMl3Pz6NV9a18NF7/sTSLa1el5V0YjHCPxvY4Jzb5JzrBR4H5g5p8wXgPudcG4BzrjkG/YqIHGBm/O/za/nVlz5EwO/jmvsX8ZM/biQS0RTPfrEI/EqgYdB+Y/TYYCcDJ5vZ62a2yMzmDHciM7vZzOrNrL6lpSUGpYlIujmtqpCnv3IBc6aX86Pfr+GWny+jq7ff67KSQiwCf7gbWA/9lRoAJgEXA9cBD5hZ4SHf5Nx851ydc66utFS3SBWRDyY/mMG9n53Jtz82ledX7uDqn7xBY5tuxBaLwG8EqgftVwFNw7T5rXOuzzm3GVjLwC8AEZG4MDP+8sITeOims2hs7WLuva9T/156z+vHIvCXAJPMrNbMMoF5wIIhbX4DXAJgZiUMTPFsikHfIiJHdPHkMn59y/nkBQNc91+LeGJJw9G/KUWNOPCdc/3ArcDzwGrgCefcSjP7vpldFW32PLDbzFYBLwPfdM7tHmnfIiLH4qSyMfz2lgs4p3Yst/3qHe5+cV1artfX/fBFJG30hyPc/tS7PLm0kRvOncD3rpqO35dan6N7pPvh6xMFRCRtBPw+/vnTpzF2TCb3v7KJ1n29/Ou1p5MV8HtdWkIo8EUkrZgZ37piKiW5Wdz57Gr2dPdy/w11jMlK/TjUvXREJC194cMncNdnTmfRplauf+BNOkJ9XpcUdwp8EUlbV59ZxU+uP4OV29q58aHFdKZ46CvwRSStfWR6Ofd+9gzebWznpp8uYW9P6l6Vq8AXkbQ355Ry/v26mSxv2MNNDy1O2dBX4IuIAFecWsG/zZvJWw17+IufLiHUF/a6pJhT4IuIRH3stAp+fO0Mlmxp5Su/eItwit1pU4EvIjLIx08fz3evnMYLq3by3QUrUuqK3NRfeCoicpxuOr+WHR09/OcrGynPD3LrrNS416MCX0RkGLddPpmdHSH+5YV1lOUHuaau+ujflOQU+CIiw/D5jB9dfRq79vbwrafepaoomw+dWOJ1WSOiOXwRkcPIDPj4j+vPoLYkl1t//tao/xAVBb6IyBHkBTOYf8OZ9IUj3PyzpXT3jt7lmgp8EZGjOKF0DP82byard3Twt796Z9Su3FHgi4gcg0umlPHNyyez4O0m5r86Oj+wT4EvInKMvnTRiXzs1Ar+6fm1LN0y+j4fV4EvInKMzIwfXn0qlYXZfOUXy2nvHl1311Tgi4gch7xgBvfMm8HOjhB3/PrdUTWfr8AXETlOM2uK+NplJ/P0O9v55dJGr8s5Zgp8EZEP4IsXnch5J4zlewtWsqllr9flHBMFvojIB+D3GXdfO4PMgI/bnnyHyCi4s2ZMAt/M5pjZWjPbYGa3H6Hdp83MmVldLPoVEfFSeUGQ71w5jfotbTz65havyzmqEQe+mfmB+4ArgGnAdWY2bZh2ecBXgDdH2qeISLL45MxKLjq5lB89tybpb70QixH+2cAG59wm51wv8Dgwd5h2PwD+CQjFoE8RkaRgZtz5yVNwwB2/Tu7758ci8CuBhkH7jdFjB5jZTKDaOff0kU5kZjebWb2Z1be0tMSgNBGR+KsqyuFv50zhlXUtPLVsm9flHFYsAt+GOXbgV5yZ+YC7ga8f7UTOufnOuTrnXF1paWkMShMRSYwbzp3AmROK+H/PrqYjlJwXZMUi8BuBwZ8MUAU0DdrPA04B/mhm7wHnAgv0xq2IpBKfz/iHq6bT2tXLvy9c73U5w4pF4C8BJplZrZllAvOABfufdM61O+dKnHMTnXMTgUXAVc65+hj0LSKSNE6pLOCaM6v56evvJeXa/BEHvnOuH7gVeB5YDTzhnFtpZt83s6tGen4RkdHkG5dPJpjh585nVntdyiFi8hGHzrlngWeHHPvOYdpeHIs+RUSSUWleFl+edRL/+Nwa/ri2mYsnl3ld0gG60lZEJMZuOn8iE8fm8I/PrkmqK3AV+CIiMZYV8PO1y05m7c5Onnl3u9flHKDAFxGJgytPG8/J48bw4z+sI5wko3wFvohIHPh9xldnn8zGln0seDs5LsZS4IuIxMmc6eVMKc/jnj+spz8c8bocBb6ISLz4fMbfXHYy7+3u4jfLm47+DfGux+sCRERS2WXTxjGlPI//enWT5zdWU+CLiMSRmfH5C2pZu7OT1zbs8rQWBb6ISJxdNWM8pXlZPPCnzZ7WocAXEYmzrICfG8+bwCvrWli3s9OzOhT4IiIJ8NlzJhDM8PGgh6N8Bb6ISAIU52byqTOq+M3ybbR3e3O/fAW+iEiCXHdWDT39ERa87c0STQW+iEiCnFKZz9SKfJ5Y0nD0xnGgwBcRSRAz45q6Kt7d1s6qpo6E96/AFxFJoE/MqCTT7+OJ+sSP8hX4IiIJVJSbyUemj+PXb22jtz+x99dR4IuIJNinzqikvbuP1zcm9spbBb6ISIKdf1IJeVkBnkvwh6Mo8EVEEiwr4Gf2tHG8sGonfQm8bbICX0TEA1ecUs6erj7e2Lg7YX0q8EVEPPDhk0vJzfTz3IrETevEJPDNbI6ZrTWzDWZ2+zDP/42ZrTKzd8xsoZlNiEW/IiKjVTDDzyVTynhx1U4iCfrM2xEHvpn5gfuAK4BpwHVmNm1Is7eAOufcacCTwD+NtF8RkdHu4sll7Nrby+odibkIKxYj/LOBDc65Tc65XuBxYO7gBs65l51zXdHdRUBVDPoVERnVLpxUAsBr6xOzPDMWgV8JDL5krDF67HA+Dzw33BNmdrOZ1ZtZfUtLSwxKExFJXuPyg0wel8efRlHg2zDHhp2QMrPPAXXAPw/3vHNuvnOuzjlXV1paGoPSRESS24WTSlj8XivdveG49xWLwG8EqgftVwGH3PvTzGYDdwBXOed6YtCviMiod8GkEnr7Iyx+rzXufcUi8JcAk8ys1swygXnAgsENzGwmcD8DYd8cgz5FRFLCObVjCfiMxZvjvx5/xIHvnOsHbgWeB1YDTzjnVprZ983sqmizfwbGAL80s+VmtuAwpxMRSSvZmX6mVuTz1tY9ce8rEIuTOOeeBZ4dcuw7g7Znx6IfEZFUNLOmkF8tbSQccfh9w70tGhu60lZExGMzawrZ1xtmfXNnXPtR4IuIeGxmdRFA3Kd1FPgiIh6bMDaH4txMlm1pi2s/CnwREY+ZGdPH57Nmh6Z0RERS3qSyPDY0743rjdQU+CIiSWDSuDF094XZtqc7bn0o8EVEksDJ48YAsG5n/KZ1FPgiIkngpLI8ANbt3Bu3PhT4IiJJoCA7g/L8YFzX4ivwRUSSRHVxNtvaNIcvIpLyKgqy2d4eitv5FfgiIklifGE229u747Y0U4EvIpIkxhcG6Qs7du2Lz0eGKPBFRJJERUE2AE174jOto8AXEUkSFQVBAHa0x+eNWwW+iEiSKMzJAKC9uy8u51fgi4gkiYJsBb6ISFoYkxXAZ9DR3R+X8yvwRUSShJmRn52hEb6ISDooyM6gI6TAFxFJeWOyAuwNaUpHRCTlBXxGfzJfaWtmc8xsrZltMLPbh3k+y8z+J/r8m2Y2MRb9ioikGr/PCCdr4JuZH7gPuAKYBlxnZtOGNPs80OacOwm4G/jRSPsVEUlFAZ8veQMfOBvY4Jzb5JzrBR4H5g5pMxd4JLr9JHCpmVkM+hYRSSk+H0kd+JVAw6D9xuixYds45/qBdmDs0BOZ2c1mVm9m9S0tLTEoTURkdAn4fIRd8gb+cCP1odUeSxucc/Odc3XOubrS0tIYlCYiMrq4Q6MxZmIR+I1A9aD9KqDpcG3MLAAUAK0x6FtEJKWE+iIEM+KzgDIWZ10CTDKzWjPLBOYBC4a0WQDcGN3+NPCSc3H6m0VEZBQL9YUJBvxxOXdgpCdwzvWb2a3A84AfeMg5t9LMvg/UO+cWAA8C/21mGxgY2c8bab8iIqko1BcmmJmkgQ/gnHsWeHbIse8M2g4Bn4lFXyIiqSzUF4nbCF9X2oqIJJFQXzip5/BFRCQGnHN0hvoZE4zJ5MshFPgiIkmiI9RPbzhC6ZisuJxfgS8ikiR27+0BYOyYzLicX4EvIpIkdu/rBaBEI3wRkdR2YISfq8AXEUlpLXsHRvia0hERSXGNbV1k+E1TOiIiqW7r7i6qi3Pw++Jz93gFvohIktiyu4sJxTlxO78CX0QkCTjn2NraxYSxuXHrQ4EvIpIEdu/rZW9PPzUa4YuIpLYtu/cBMGGsAl9EJKWt2t4JwJSK/Lj1ocAXEUkCq5raKczJYHxBMG59KPBFRJLAqqYOplXkYxafJZmgwBcR8Vx/OMKaHZ1MHx+/6RxQ4IuIeG7Trn309EeYpsAXEUltbzfsAeCU8QVx7UeBLyLiscWbWynKyeDE0jFx7UeBLyLiscXvtXLWxGJ8cbqHzn4KfBERD+1oD7Fldxdn1xbHva8RBb6ZFZvZi2a2PvpYNEybGWb2hpmtNLN3zOzakfQpIpJKFr/XCsA5tWPj3tdIR/i3Awudc5OAhdH9obqA/+Wcmw7MAX5sZoUj7FdEJCUs3ryb3Ew/Uyvy4t7XSAN/LvBIdPsR4BNDGzjn1jnn1ke3m4BmoHSE/YqIjHrOOV5Z18I5J4wl4I//DPtIexjnnNsOEH0sO1JjMzsbyAQ2Hub5m82s3szqW1paRliaiEhy29iyl4bWbi6desTojJnA0RqY2R+A8mGeuuN4OjKzCuC/gRudc5Hh2jjn5gPzAerq6tzxnF9EZLRZuLoZgFlTkiTwnXOzD/ecme00swrn3PZooDcfpl0+8Azwbefcog9crYhIClm4pplpFflUFGQnpL+RTuksAG6Mbt8I/HZoAzPLBH4N/Mw598sR9icikhL2dPWydEtbwqZzYOSB/0PgMjNbD1wW3cfM6szsgWiba4APAzeZ2fLo14wR9isiMqq9vLaZcMRxSYKmc+AYpnSOxDm3G7h0mOP1wF9Gtx8FHh1JPyIiqea3y5uoLMxmRlXiVqnrSlsRkQTbvbeHP63fxVUzxsf9dgqDKfBFRBLsmXe3E4445s4Yn9B+FfgiIgn22+VNTB6Xx5Ty+N7/figFvohIAjW0drF0SxtzZyZ2dA8KfBGRhHqivgEzmDujMuF9K/BFRBKktz/CLxY3MGtyGZWFibnYajAFvohIgjy/cge79vbwufMmeNK/Al9EJEEeXbSF6uJsLprkzQ2DFfgiIgmwbmcnb25u5XPnTEjo2vvBFPgiIgnw8J/fIzPg4zN11Z7VoMAXEYmz5o4QT9Y3cvUZVRTnZnpWhwJfRCTOHnxtM/2RCF+86ARP61Dgi4jEUXtXH48u2sKVp41nwthcT2tR4IuIxNEjb7zHvt4wX7r4RK9LUeCLiMTLvp5+fvr6ZmZNKWNqRWLvmzMcBb6ISJw88KfNtHX18eVZJ3ldCqDAFxGJi117e5j/6kauOKWcmTVFXpcDKPBFROLi3pc2EOqP8I3LJ3tdygEKfBGRGNuyex+PvbmFa8+q5sTSMV6Xc4ACX0Qkxv7lhXUEfD6+eukkr0s5iAJfRCSGFm3aze/ebuILF9ZSlh/0upyDKPBFRGKkLxzh73+zgqqibL50cXKszBlsRIFvZsVm9qKZrY8+HvataDPLN7NtZnbvSPoUEUlWD722mfXNe/nex6eTnen3upxDjHSEfzuw0Dk3CVgY3T+cHwCvjLA/EZGk1LSnm3sWrmf21HHMnjbO63KGNdLAnws8Et1+BPjEcI3M7ExgHPDCCPsTEUlKP3h6FRHn+O7Hp3ldymGNNPDHOee2A0Qfy4Y2MDMfcBfwzaOdzMxuNrN6M6tvaWkZYWkiIonx9DtNPLdiB1+eNYnq4hyvyzmswNEamNkfgPJhnrrjGPv4K+BZ51yD2ZE/5cU5Nx+YD1BXV+eO8fwiIp5p7gzx979ZwelVBfyfD3t7++OjOWrgO+dmH+45M9tpZhXOue1mVgE0D9PsPOBCM/srYAyQaWZ7nXNHmu8XEUl6zjn+7qkVdPWGueuaGQT8yb3wcaTVLQBujG7fCPx2aAPn3PXOuRrn3ETgG8DPFPYikgqeWraNP6zeyTcvn8xJZclzRe3hjDTwfwhcZmbrgcui+5hZnZk9MNLiRESS1bY93Xzvdys5u7aYvzi/1utyjslRp3SOxDm3G7h0mOP1wF8Oc/xh4OGR9Cki4rXe/gi3/nwZzsG/fPp0fL4jvz+ZLEYU+CIi6ehHv1/DW1v38B/Xn0HN2ORdlTNUcr/DICKSZH6/YgcPvraZmz40kY+eWuF1OcdFgS8icoy27N7HN598m9OrC/m7j071upzjpsAXETkGob4wt/x8GT4z7r1uJpmB0RefmsMXETkGzsHJZXl89dKTk/pq2iNR4IuIHIPsTD//eu0Mr8sYkdH3N4mIiHwgCnwRkTShwBcRSRMKfBGRNKHAFxFJEwp8EZE0ocAXEUkTCnwRkTRhziXnJwmaWQuw5Ti+pQTYFadyko1ea+pJl9cJeq3xNsE5VzrcE0kb+MfLzOqdc3Ve15EIeq2pJ11eJ+i1eklTOiIiaUKBLyKSJlIp8Od7XUAC6bWmnnR5naDX6pmUmcMXEZEjS6URvoiIHIECX0QkTaRU4JvZZ8xspZlFzCxplkLFipnNMbO1ZrbBzG73up54MrOHzKzZzFZ4XUs8mVm1mb1sZqujP7t/7XVN8WJmQTNbbGZvR1/rP3hdUzyZmd/M3jKzp72uZb+UCnxgBfAp4FWvC4k1M/MD9wFXANOA68xsmrdVxdXDwByvi0iAfuDrzrmpwLnALSn879oDzHLOnQ7MAOaY2bke1xRPfw2s9rqIwVIq8J1zq51za72uI07OBjY45zY553qBx4G5HtcUN865V4FWr+uIN+fcdufcsuh2JwMBUeltVfHhBuyN7mZEv1Jy1YiZVQEfAx7wupbBUirwU1wl0DBov5EUDYZ0ZWYTgZnAm95WEj/RaY7lQDPwonMuVV/rj4HbgIjXhQw26gLfzP5gZiuG+UrZ0W6UDXMsJUdH6cjMxgC/Ar7qnOvwup54cc6FnXMzgCrgbDM7xeuaYs3MrgSanXNLva5lqIDXBRwv59xsr2vwSCNQPWi/CmjyqBaJITPLYCDsH3POPeV1PYngnNtjZn9k4H2aVHtj/nzgKjP7KBAE8s3sUefc5zyua/SN8NPYEmCSmdWaWSYwDws2x+oAAAC/SURBVFjgcU0yQmZmwIPAaufcv3pdTzyZWamZFUa3s4HZwBpvq4o959y3nHNVzrmJDPx/+lIyhD2kWOCb2SfNrBE4D3jGzJ73uqZYcc71A7cCzzPwxt4TzrmV3lYVP2b2C+ANYLKZNZrZ572uKU7OB24AZpnZ8ujXR70uKk4qgJfN7B0GBjAvOueSZsliOtCtFURE0kRKjfBFROTwFPgiImlCgS8ikiYU+CIiaUKBLyKSJhT4IiJpQoEvIpIm/j/TDGPVCOzEUAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x,v = pyaudi_integration.y.T\n",
    "plt.plot([it for it in x.astype(np.float64)],[it for it in v.astype(np.float64)])\n",
    "plt.plot(x[0].constant_cf, v[0].constant_cf, 'ro')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "xf, vf = x[-1], v[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "initial xf: 3.8065257635631173\n",
      "initial vf: -0.4637962920728024\n"
     ]
    }
   ],
   "source": [
    "print(\"initial xf: {}\".format(xf.constant_cf))\n",
    "print(\"initial vf: {}\".format(vf.constant_cf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Studying the effects of the weights on the behavior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have represented all the robot behavior (x, v) as a polynomial function of the weights. So we now know what happens to the behaviour if we change the weights!! Lets see … we only consider the final state, but the same can be done for all states before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "new xf: 11.843378895251684\n",
      "new vf: 3.5577687770021136\n"
     ]
    }
   ],
   "source": [
    "dweights = dict({'da': -0.002, 'db': 0.003, 'dc': -0.02, 'de': 0.03, 'df': 0.02, 'dg': -0.01})\n",
    "#Lets predict the new final position of our 'robot' if we change his controller as defined above\n",
    "print(\"new xf: {}\".format(xf.evaluate(dweights)))\n",
    "print(\"new vf: {}\".format(vf.evaluate(dweights)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check that we learned the correct map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now simulate again our behavior using the new weights to see where we end up to check if the prediction made after our differential learning is correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_weights = np.array([it + dweights['d' + it.symbol_set[0]] for it in weights])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "79dbf937869b407e8665ba69badc068d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=9000000000.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "pyaudi_integration2 = de.OdeSystem(rhs, y0=y0, dense_output=True, t=(0, 10.), dt=0.01, rtol=1e-12, atol=1e-12, constants=dict(weights=new_weights))\n",
    "\n",
    "pyaudi_integration2.set_method(\"RK45\")\n",
    "pyaudi_integration2.integrate(eta=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdoAAAHSCAYAAABLtwrCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3hU1cLF4d9JryT0qgRQWqghIEov0kHpVUAvYr2ICoINlGpBUFAEVFSk9yIiTbkUkSZFqiiEXgMJ6clkzvdHMJ9IgASSnExmvc+Th2TmeGZNDFnsU/Y2TNNEREREsoaL1QFERERyMxWtiIhIFlLRioiIZCEVrYiISBZS0YqIiGQhFa2IiEgWcsuKnRYoUMAMCgrKil2LiIjkOLt27bpsmmbBtJ7LkqINCgpi586dWbFrERGRHMcwjBO3ek6HjkVERLKQilZERCQLqWhFRESyUJaco01LUlISp0+fJj4+PrteUhyAl5cXJUqUwN3d3eooIiJZItuK9vTp0/j7+xMUFIRhGNn1spKDmaZJeHg4p0+fplSpUlbHERHJEtl26Dg+Pp78+fOrZCWVYRjkz59fRzlEJFfL1nO0Kln5N/1MiEhup4uh/qVVq1ZERETcdpthw4axbt26u9r/hg0baNOmzV39tyIi4niy7RxtTmeaJqZp8sMPP9xx2xEjRmRDIhERyQ2cakQ7fvx4KlWqRKVKlfj4448JCwujQoUKPP/884SEhHDq1CmCgoK4fPkyACNHjqR8+fI8+uijdO/enXHjxgHQt29fFi5cCKTMgjV8+HBCQkKoXLkyhw8fBmD79u088sgjVK9enUceeYQjR45Y86ZFRMRSloxo311xgINnr2XqPisWy8PwtsG3fH7Xrl18/fXXbNu2DdM0eeihh2jQoAFHjhzh66+/ZvLkyTdsv3PnThYtWsTu3bux2WyEhIRQo0aNNPddoEABfvvtNyZPnsy4ceP48ssvKV++PBs3bsTNzY1169bxxhtvsGjRokx9zyIikvM5zaHjzZs30759e3x9fQHo0KEDmzZtomTJktSuXTvN7R977DG8vb0BaNu27S333aFDBwBq1KjB4sWLAYiMjKRPnz4cPXoUwzBISkrK7LckIiIOwJKivd3IM6uYppnm438Xb3q3T4unpycArq6u2Gw2AN5++20aNWrEkiVLCAsLo2HDhhkLLCIiuYLTnKOtX78+S5cuJTY2lpiYGJYsWUK9evVuuX3dunVZsWIF8fHxREdHs3Llygy9XmRkJMWLFwfgm2++uZfoIiLiwJzm0HFISAh9+/alVq1aAPTr14+8efPecvuaNWvSrl07qlatSsmSJQkNDSUgICDdr/faa6/Rp08fxo8fT+PGje85v4iIOCYjI4dI0ys0NNT893q0hw4dokKFCpn+WlkpOjoaPz8/YmNjqV+/PtOmTSMkJMTqWLmOI/5siIj8k2EYu0zTDE3rOacZ0d6N/v37c/DgQeLj4+nTp49KVkQkN/jzTyhTBrJpZjoV7W3Mnj3b6ggiIpKZ/vgDqleHYcNgyJBseUmnuRhKREScXFIS9OwJXl7Qq1e2vaxGtCIi4hzeeQd27oSFC+H6XSHZQSNaERHJ/TZuhLFj4amnoGPHbH1pFa2IiORuV6+mHCouUwY++STbX96pi7Zfv34cPHgwU/b1z8UIbmXMmDE3fP3II49kymuLiMgtmCY8+yycOwezZ4OfX7ZHcOqi/fLLL6lYsWK2vd6/i/aXX37JttcWEXFKM2bA/PkwYgTUrGlJBKcp2piYGFq3bk3VqlWpVKkS8+bNo2HDhvw9sYafnx9DhgyhRo0aNG3alO3bt9OwYUNKly7N8uXLgZSpFF988cXUfbZp04YNGzbc9FqPP/44NWrUIDg4mGnTpgEwdOhQ4uLiqFatGj179kx9TUiZV3nw4MFUqlSJypUrM2/ePCBlkfiGDRvSqVMnypcvT8+ePTM0B7OIiFP76y948UVo0ABee82yGNZcdTxwIOzZk7n7rFYNPv74lk//+OOPFCtWLHXO4sjISD7//PPU52NiYmjYsCHvv/8+7du356233mLt2rUcPHiQPn360K5du3RHmT59Ovny5SMuLo6aNWvSsWNH3nvvPT799FP2pPG+Fy9ezJ49e9i7dy+XL1+mZs2a1K9fH4Ddu3dz4MABihUrRp06ddiyZQt169ZNdxYREaf09608bm7w3Xfg6mpZFKcZ0VauXJl169YxZMgQNm3adNO8xR4eHrRo0SJ12wYNGuDu7k7lypUJCwvL0GtNnDiRqlWrUrt2bU6dOsXRo0dvu/3mzZvp3r07rq6uFC5cmAYNGrBjxw4AatWqRYkSJXBxcaFatWoZziIi4pRGjoRt22DqVLjvPkujWDOivc3IM6uULVuWXbt28cMPP/D666/TrFmzG553d3fHuD4dl4uLS+rSdy4uLqlL37m5uWG321P/m/j4+JteZ8OGDaxbt46tW7fi4+NDw4YN09zun253OPjvHHDjMnwiInILmzbB6NHQty906WJ1GucZ0Z49exYfHx969erFoEGD+O233zK8j6CgIPbs2YPdbufUqVNs3779pm0iIyPJmzcvPj4+HD58mF9//TX1OXd39zQXgK9fvz7z5s0jOTmZS5cusXHjxtRVhkREJAMiIuCJJ6BUKZg40eo0gBPNDPX7778zePBgXFxccHd35/PPP2fQoEEZ2kedOnUoVaoUlStXplKlSmkuMtCiRQumTJlClSpVKFeuHLVr1059rn///lSpUoWQkBBmzZqV+nj79u3ZunUrVatWxTAMPvjgA4oUKcLhw4fv/g2LiDgb04RnnoEzZ2DzZvD3tzoRoGXyJAfQz4aIZIrp0+E//4ExY+D117P1pW+3TJ7THDoWEZFc7PBh+O9/oXFjS2/lSYuKVkREHFtCAnTvDt7elt/KkxanOUcrIiK51Ouvp8zNsHw5FCtmdZqbaEQrIiKOa9UqmDABXngB2ra1Ok2aVLQiIuKYzp+HPn2gcmX48MPbbpqUnERUQlQ2BbuRilZERByP3Z5SslFRMHduyvnZWzBNk34r+lH/m/rE224/gVBWcNqifeeddxg3bhwAw4YNY926dQBs2rSJ4OBgqlWrRlxcHIMHDyY4OJjBgwdnaZ6PP/6Y2NjY1K9btWpFRETEbf+bWy3Nl54l++7Vzp07GTBgQJa+hojILU2YAGvWpPx5h1XY3vzpTWbsnUHHCh3xcvPKpoD/TxdDASNGjEj9fNasWQwaNIgnn3wSgKlTp3Lp0qUbpkK8HZvNhptbxr+tH3/8Mb169cLHxweAH374IcP7yGzJycm43uLqvdDQUEJD07xlTEQka+3YAUOHQvv2KRNU3Man2z9l7OaxPFvjWd6s92Y2BbyRU41oR48eTbly5WjatClHjhxJfbxv374sXLiQL7/8kvnz5zNixAh69uxJu3btiImJ4aGHHmLevHlcunSJjh07UrNmTWrWrMmWLVuAlNFx//79adasGb179yY5OZnBgwdTs2ZNqlSpwtSpU4FbL3s3ceJEzp49S6NGjWjUqBFw46g0rWX30mvmzJnUqlWLatWq8cwzz5CcnAzAc889R2hoKMHBwQwfPjx1+6CgIEaMGEHdunVZsGABDRs2ZMiQIdSqVYuyZcuyadOm1PfSpk2b1Pf/1FNPpS4rOPEf056NHDmS8uXL8+ijj9K9e/fUowgiInclMhK6dUu5uvirr+D6HPVpWXxoMQNWDeCxco/xaatPU+ezz26WjGgH/jiQPeczd5m8akWq8XGLWy9WsGvXLubOncvu3bux2WyEhIRQo0aNG7bp168fmzdvpk2bNnTq1AlIWTP276XtevTowcsvv0zdunU5efIkzZs359ChQ6n737x5M97e3kybNo2AgAB27NhBQkICderUSV3EIK1l7wYMGMD48eP5+eefKVCgwE3Z01p2L3/+/Hf8nhw6dIh58+axZcsW3N3def7555k1axa9e/dm9OjR5MuXj+TkZJo0acK+ffuoUqUKAF5eXmzevBmAKVOmYLPZ2L59Oz/88APvvvtu6mH2fzp8+DA///wzUVFRlCtXjueee469e/eyaNGi237PRUTS7e8pFk+cgI0bIW/eW2666cQmeizqwcP3PcycjnNwdbHu3lqnOXS8adMm2rdvn3poNiPry/5t3bp1HDx4MPXra9euERUVlbo/7+sn49esWcO+fftYuHAhkLLQwNGjR/Hw8Ehd9g5IXfbuTuvLTpw4kSVLlgCkLruXnqJdv349u3btombNmgDExcVRqFAhAObPn8+0adOw2WycO3eOgwcPphZt165db9hPhw4dAKhRo8Ytl+lr3bo1np6eeHp6UqhQIS5cuMDmzZt57LHHUr8vbXPopfci4iC++grmzUuZYvGRR2652YGLB2g3tx1BgUEs77Ycb/dbXyiVHSwp2tuNPLPSvR42sNvtbN26NbU4/snX1zf1c9M0mTRpEs2bN79hmw0bNmR42bu7WXbvnzn69OnD2LFjb3j8+PHjjBs3jh07dpA3b1769u17wz7/+V7g/5fqu13etN5XVsyjLSJO6sABGDAAmjaFIUNuudnpa6dpMasFXm5e/NjrR/L73HlQktWc5hxt/fr1WbJkCXFxcURFRbFixYoM76NZs2Z8+umnqV//fUj535o3b87nn3+euiTeH3/8QUxMzG337e/vnzo6/qfbLbt3J02aNGHhwoVcvHgRgCtXrnDixAmuXbuGr68vAQEBXLhwgVWrVqV7nxlRt25dVqxYQXx8PNHR0axcuTJLXkdEcrnYWOjaNWU1nu++A5e0qysiPoKWs1oSGR/Jqp6rCAoMyt6ct+A0h45DQkLo2rUr1apVo2TJktSrVy/D+5g4cSIvvPACVapUwWazUb9+faZMmXLTdv369SMsLIyQkBBM06RgwYIsXbr0tvvu378/LVu2pGjRovz888+pj99u2b07qVixIqNGjaJZs2bY7Xbc3d357LPPqF27NtWrVyc4OJjSpUtTp06d9H8TMqBmzZq0a9eOqlWrUrJkSUJDQwkICMiS1xKRXOzll1NGtKtXQ5EiaW6SYEvg8bmPc+TyEVb1XEW1ItWyOeStaZk8yVLR0dH4+fkRGxtL/fr1mTZt2k3r+OpnQ0Ruaf78lNHskCHw3ntpbmI37XRf1J35B+Yzq8MselTukc0hb79MntOMaMUa/fv35+DBg8THx9OnT5+bSlZE5JaOHYOnn4batWHkyDQ3MU2TV1a/wvwD8/nw0Q8tKdk7UdFKlpo9e7bVEUTEESUkpIxkXVxgzhxwd09zs4+2fsQn2z5h4EMDefXhV7M5ZPqoaEVEJOcZMgR27oQlSyAoKM1NZv8+m8FrB9MluAsfNf/Isgkp7iRbrzrW7R7yb/qZEJGbLFkCn3wCL70Ejz+e5ibrj62n79K+NAxqyIzHZ+Bi5NybaLItmZeXF+Hh4frFKqlM0yQ8PBwvr+yf5FtEcqjjx+GppyA0FD74IM1N9pzfQ/t57SlfoDxLui7B0y19c9FbJdsOHZcoUYLTp09z6dKl7HpJcQBeXl6pM2WJiJNLTEw5L2uaKTNAeXjctElYRBgtZ7Uk0CuQVT1XEegVaEHQjMm2onV3d6dUqVLZ9XIiIuJohg5NWZln0SIoXfqmp8Njw2kxswXxtnjWPbmO4nmKWxAy43QxlIiIWG/ZspS1ZV98Ea7Pr/5PsUmxtJ3TlrCIMNY+sZbgQsEWhLw7KloREbFWWBj07QshIZDGUpo2u43ui7rz6+lfWdhlIfVKZnxmPyupaEVExDqJiSnry9rtKbNAed54YZNpmryw8gWWH1nOpJaT6FDh5tFuTqeiFRER67z2GmzbBgsWQJkyNz09auMopv02jaF1hvJirRctCHjvcu6NRyIikrstWpRyv+yAAdCp001PT989nWEbhtG7am/GNBljQcDMoaIVEZHs9+efKffL1qoFH35409Mr/1hJ/xX9aVamGV+2/TLHzvqUHipaERHJXnFxKSNYV9eU87L/ul92+5ntdFnYhapFqrKw80LcXdOe59hR6BytiIhkr5degr174fvvoWTJG546Gn6U1rNbU9i3MCt7rMTf09+ikJlHI1oREck+330HX3yRMjlF69Y3PHUh+gItZrUAYHWv1RTxS3uRd0ejEa2IiGSPAwfg2WehQYOb1peNToym9ezWnI8+z0+9f+LB/A9aFDLzqWhFRCTrRUdD587g55eyvqzb/9dPUnISneZ3Ys/5PSzrtoyHSjxkYdDMp6IVEZGsZZrw9NNw5AisXQtFi/7jKZN+K/qx+q/VfNH2C1qXbX2bHTkmFa2IiGStzz6DuXNhzBho3PiGp9766S1m7J3Buw3fpV9IP4sCZi1dDCUiIlln61Z45RVo0waGDLnhqck7JjNm8xieDnmat+u/bVHArKeiFRGRrHHpEnTpAiVKwIwZ4PL/lbP40GJe/OFF2pVrx+TWkx16Qoo70aFjERHJfMnJ0KNHStlu3Qp586Y+tfnkZnos6sFDJR5iTsc5uLnk7irK3e9ORESs8c47sG4dfPklVK+e+vDBSwdpO6ctJQNLsqL7CnzcfazLmE106FhERDLXypUwalTKXMb/+U/qw2eunaHFzBZ4uXmxutdqCvgUsDBk9tGIVkREMs/x4/DEE1CtGnz6aerDEfERtJzVkoj4CP7X938EBQZZlzGbqWhFRCRzxMVBx44pi7gvXAje3gAk2BJoP689hy4fYlXPVVQvWv0OO8pdVLQiInLvTBOefx5274YVK1IXcbebdnov7c2GsA3MbD+TpqWbWhw0++kcrYiI3LupU+Gbb2DYsJR7Zq8btGYQ8w/M54OmH9CzSk/r8llIRSsiIvfm119hwABo2RKGD099+KNfPmLCrxMYUGsAgx4ZZGFAa6loRUTk7l24kHJetkQJmDkzdVKKOb/PYdDaQXSu2JkJLSbk6gkp7kTnaEVE5O7YbNC1K1y5kjIpRb58AKw/tp4+S/vQoGQDZrSfgYvh3GM6Fa2IiNydIUPgf/9LmV6xWjUA9p7fS/t57SmbvyxLuy3Fy83L4pDWc+5/ZoiIyN2ZPx/Gj4cXXki5bxYIiwij5ayWBHgF8GOvHwn0CrQ4ZM6gEa2IiGTM77/Dk0/Cww+nlC0QHhtOi5ktiLPFsfnJzZTIU8LikDmHilZERNLvyhV4/HEICIBFi8DDg7ikONrOaUtYRBhrnlhDcKFgq1PmKCpaERFJn79X5Dl1KuXcbNGi2Ow2ui/qzq+nf2V+5/nUL1nf6pQ5jopWRETS5+23YfXqlMkpHn4Y0zR58YcXWXZkGRNbTKRTxU5WJ8yRdDGUiIjc2cKFMHYsPP009O8PwKiNo5i6aypD6gzhvw/91+KAOdcdi9YwDC/DMLYbhrHXMIwDhmG8mx3BREQkh9i/H/r2hdq1YdIkAL7Y9QXDNgzjiSpPMKbJGGvz5XDpOXScADQ2TTPaMAx3YLNhGKtM0/w1i7OJiIjVrl5NufjJ3z/l4idPT5YdXsazK5+lxQMt+KrdV04/IcWd3LFoTdM0gejrX7pf/zCzMpSIiOQAycnQsyecPAk//wzFirH55Ga6LepGaLFQFnRegLuru9Upc7x0/TPEMAxXwzD2ABeBtaZpbsvaWCIiYrm33oJVq1IOF9epw/6L+2k7py33B9zPyh4r8fPwszqhQ0hX0ZqmmWyaZjWgBFDLMIxK/97GMIz+hmHsNAxj56VLlzI7p4iIZKe5c+G99+CZZ+CZZzgZeZIWM1vg7ebN6l6rKeBTwOqEDiNDB9ZN04wANgAt0nhummmaoaZphhYsWDCT4omISLbbvRueegrq1oWJE1NnfYpKjOLHXj8SFBhkdUKHkp6rjgsahhF4/XNvoClwOKuDiYiIBS5eTLn4KX9+WLiQWMNGmzltOHb1GMu7LadK4SpWJ3Q46bnquCjwrWEYrqQU83zTNL/P2lgiIpLtkpKgc+eUst28GVvB/HSd155tp7exoPMCGgQ1sDqhQ0rPVcf7gOrZkEVERKw0cCBs3AizZmGGhNB/+X/4/o/vmdxqMh0rdrQ6ncPSzU8iIgJffAGTJ8PgwdCjB2/+9CZf7/maYfWH8VzN56xO59BUtCIizm7TppR1ZZs3h7FjmbRtEmM3j6V/SH/eafiO1ekcnopWRMSZhYVBhw5QujTMncv8w4t46ceXeKzcY3zW+jMMw7A6ocNT0YqIOKvoaHjsMbDZYPly1l/ZRa/Fvahzfx3mdJyDm4sWeMsM+i6KiDgjux2eeCJlwYBVq9jtH0P7b9pTrkA5lndbjre7t9UJcw0VrYiIMxo+HJYuhY8/5ljNB2j51SMEegWyqucq8nrntTpdrqKiFRFxNnPnwqhR8J//cOGprjT7ui5J9iR+7vMzJfKUsDpdrqOiFRFxJrt2wZNPQt26RE14n1ZzmnE26iw/9fmJCgUrWJ0uV1LRiog4izNnoF07KFSIxPlz6LC0G3vP72VZt2XULlHb6nS5lopWRMQZxMamXGF87Rr2zZvo+8trrDu2juntptO6bGur0+VqKloRkdzOboc+feC33zCXLmXQ+RnM2T+HsU3G8mT1J61Ol+vpPloRkdzunXdg4UL48EPez3eQCb9OYECtAQypM8TqZE5BI1oRkdxs9mwYORKeeoovGvrz+vfP0L1Sdya0mKBZn7KJilZEJLfaujVlAff69Vn0UjOeXdqDFg+04JvHv8HF0AHN7KKiFRHJjU6eTFnAvXhxfvpkID2Wd+Oh4g+xsPNCPFw9rE7nVFS0IiK5zbVr0KYNxMezc9EkHlvVm7L5y/J9j+/x9fC1Op3TUdGKiOQmNht06wYHD3J40VRabn2BAj4FWN1rNfm881mdzimpaEVEcpOXX4ZVqzg1eSzNjr2Li+HCml5rKOZfzOpkTktFKyKSW0yaBJ9+Svirz9PcPoPIhEg29NnAg/kftDqZU1PRiojkBitXwsCBRLdvTatyOzh24Rire62metHqVidzeipaERFHt3cvdOtGQvXKdGgXx86Tu1jcZTENghpYnUzQzFAiIo7t7Flo04bkwDz0HnA/a0/8xJdtv+Sx8o9ZnUyuU9GKiDiq6Gho0wbz6hVeGPkw84+v4MNHP9T8xTmMilZExBH9fRvP3r0MG9+OqScWMaTOEAY9MsjqZPIvKloREUdjmjBwIKxcyScfdmTUubn8p/p/GNtkrNXJJA0qWhERR/Pxx/DZZ8wc2oqBUQtoX749U9pM0SIBOZSuOhYRcSRLlsCrr7Ky7yM86bOGRvc3YnbH2bi56Nd5TqX/MyIijmLbNujZky3NK9L5gd1UKViFpd2W4uXmZXUyuQ0VrYiII/jrL2jbln0V8tGmwWnu87uPVT1Xkcczj9XJ5A5UtCIiOd3ly9CyJcf8kmjeBXzd/VjTaw2FfAtZnUzSQRdDiYjkZLGx0LYtZ66epOmzPiQayazutZqSgSWtTibppKIVEcmpkpOhZ08u7/uVR18tyKXka/zY80eCCwVbnUwyQIeORURyouv3ykauWkrzt0pwPPkyP/b8kZrFa1qdTDJIRSsikhN99BGxUz+lzZDi7DPPs6zrMi0S4KBUtCIiOc2cOSQMHUyHl4vwi/s55nSYQ6sHW1mdSu6SilZEJCdZvx5b3970fKYAq/3O81Xbr+gS3MXqVHIPdDGUiEhOsWcP9vaP06+HH4sKXWZC8wk8Vf0pq1PJPVLRiojkBGFhmC1bMLAFfBsUwbsN32Vg7YFWp5JMoEPHIiJWCw+HFi14u3oEk4ITeKX2K7xd/22rU0kmUdGKiFgpNhbatOGDwn8y+qFk+lXvx7hm47QSTy6iQ8ciIlZJSoKuXZli+5UhjZPpGtxVy93lQipaERErmCb078/sE9/zfBuD1g+25rv23+Hq4mp1MslkKloRESu8/jrLtn5D7w4GDYIasKDzAtxd3a1OJVlARSsikt0mTGD9/Pfp0s2FGsVrsrzbcrzdva1OJVlERSsikp1mzmTr+Fd4rJcrZQtVZFWvVfh7+ludSrKQilZEJLv8+CN7h/SlVR83ihYIYm3vteTzzmd1KsliKloRkeywdStHnm5Ps94GfnkLs673eor4FbE6lWQDFa2ISFbbv59j3VvQpHsS5A1kXZ/1WrjdiWjCChGRrHT8OCc7NKFxpxji8vmzoc9PlCtQzupUko1UtCIiWeX8ec62a0ST1peJyOfL+t7rqFy4stWpJJupaEVEskJEBBfbNaFJo5Ocz+/N2t5rqFGshtWpxAIqWhGRzBYbS3iHFjQNPciJgp78+MQqapeobXUqsYiKVkQkMyUmEtGlHc0f3MYfRdz5vuf31C9Z3+pUYiEVrYhIZklOJqp3V1oVWc++Yq4s6baEpqWbWp1KLKaiFRHJDHY7sc88RRvvpWy/z4UFXRbQumxrq1NJDqD7aEVE7pVpEv/qSzyeOIPNJQ1mdpxF+wrtrU4lOYSKVkTkHiWOGEani5+ytgxMf3w63Sp1szqS5CAqWhGRe5A0fhzdDo1iZVmY0moyfar1tTqS5DAqWhGRu5Q0ZTLdtw1mSQWY2GwCz9R8zupIkgOpaEVE7oJtxjf0/OkFFlWECU0+5L8PD7Q6kuRQKloRkQyyLZzPEyueZEEwjGs0loF1B1kdSXIwFa2ISAYkr1xBn3ndmFsJ3q83glfrD7U6kuRwKloRkXRKXr+WJ79+nNmVTMbWGc5rjd+2OpI4ABWtiEg62Df+j36ft+S7ynZG1X6DoU3fsTqSOAgVrYjIHdh/2cLTE5vyTeVk3q35Gm82H211JHEgKloRkduw79jOsx81YnplG8NCXmZYq/etjiQORkUrInIL9j27eeG9enxRJYk3q/6Xd9p8ZHUkcUBaVEBEJA32/b/z/OhHmFolkaGVnmPkY59gGIbVscQBqWhFRP7FfmA/z4ysxZeV4nkj+DlGdfhMJSt3TYeORUT+IfngAfqNrMmXFeN5O/h5RnVUycq90YhWROS65EMHeWpUKDMqxPNO8IsM7zTJ6kiSC6hoRUSA5COH6TsqlJnl4hkR/CJvq2Qlk6hoRcTp2f44TO9RIcwpG8fo4AG80ekTqyNJLqKiFRw6EJEAACAASURBVBGnlnTkIL1GhzL/gTjeC36JIZ0+tjqS5DIqWhFxWklHDtJ9TA0WlYlnXOVXebXDOKsjSS6kohURp5R4aD/d3q/JktLxjK8ymJfbf2B1JMmlVLQi4nTiD/1O5w9r8n2pBD6pOpQBj4+1OpLkYipaEXEqMft/4/FPHmZdyUSmhAzjmbbvWh1JcjkVrYg4jajd22g9pR5biiXxTeho+rR+w+pI4gRUtCLiFCJ2bKLF9MbsLGJjdt3xdH30ZasjiZNQ0YpIrnd563qazWzOgYLJLGr4OY81etbqSOJEVLQikqud37SKpvPb8Fc+k2XNvqZF3b5WRxIno6IVkVzr9LrFNFnRiTMB8EPrWTSq3d3qSOKEVLQikisdXzmTJut6E+5nsLr9QurUaG91JHFSKloRyXUOL5rKo788S4yPK+u7rSS0cnOrI4kTu+N6tIZh3GcYxs+GYRwyDOOAYRgvZUcwEZG7sXvWOOpve5YkDzc29P5JJSuWS8+I1ga8aprmb4Zh+AO7DMNYa5rmwSzOJiKSIZu/Gk7rP0cQ6OLBuv6/8GCpGlZHErlz0ZqmeQ44d/3zKMMwDgHFARWtiOQYP04cQIeLk7jf7s3al3ZxX7EKVkcSATJ4jtYwjCCgOrAtK8KIiNyNBe89Qc/YmQQn5mH1a/soVKCk1ZFEUt3xHO3fDMPwAxYBA03TvJbG8/0Nw9hpGMbOS5cuZWZGEZG0mSbTh7ejW9xMaiUU4Oe3jqpkJcdJV9EahuFOSsnOMk1zcVrbmKY5zTTNUNM0QwsWLJiZGUVEbma3M+H1hvzHZQWPJhRnzbvHCMxTyOpUIjdJz1XHBvAVcMg0zfFZH0lE5PbMpCSGvxrCK94b6ZT4AMtH/omPl7/VsUTSlJ4RbR3gCaCxYRh7rn+0yuJcIiJpSo6N4b8vl2NE4F6eojpzRxzCw93L6lgit5Seq443A0Y2ZBERua2Eq5fp/VYw8wtd5FWvxnz42jpSDrqJ5FzpvhhKRMRKUefCaPN2GeYXusgH+boxbsh6law4BE3BKCI53qU/9tDqs4fZnT+er0u+RN++H1sdSSTdVLQikqOd2LmOZrNacjKPjSVVxtC24+tWRxLJEBWtiORY+9fNpvnqXsR6w9qG06nb5EmrI4lkmIpWRHKkLQvH02bXq/i4uLKx4woq12hpdSSRu6KiFZEcZ/m0V+h2cgIlbJ6s6b+JoAdrWh1J5K6paEUkR/n8vY68GLeYGnH+fD/oNwoVfcDqSCL3REUrIjmCPdnGm289wnteO2gTU5S57+zH1z+f1bFE7pmKVkQslxBzjafersLsgBM8G1+ZSe/txM3Nw+pYIplCE1aIiKUiLpyg5VslmR1wgrFuLZg8eo9KVnIVjWhFxDKnDm2j1RcNOOKfwHcFn6HX81OsjiSS6VS0ImKJfZsX0WpZF6K87Kyq+gFNOg62OpJIllDRiki2W7foQzrseo08uLC55UIq1+todSSRLKOiFZFs9d1nz/DUhWmUj/dkVf+NlChfy+pIIllKRSsi2cK02xkzoilvGT/TKCovS4buJqBwSatjiWQ5Fa2IZLmkuBheeKsaX+T5k57RpZk+Zi8e3n5WxxLJFrq9R0Sy1NVzx2nxegm+yPMnb1Cf794/qpIVp6IRrYhkmaO719NmZkuO50ni20LP0Ps53b4jzkdFKyJZYsOKSXT45SVc3GF96ETqtfuv1ZFELKGiFZFM99Xkfjx7/iseTPBgxRM/UKZ6E6sjiVhGRSsimSbZlsTrIxvwoctWmkXmZf7Q3wgoGmR1LBFLqWhFJFNER1yk18jqLMtzluejyvPJ+7tw8/KxOpaI5VS0InLPTv+xg7ZTGrDPP45JHo/z4oeLwTCsjiWSI6hoReSe7Fj3LY+teYpoLzsry75Di57DrY4kkqOoaEXkrs2f9hJ9Tk6kiN2Vta0WEVy3g9WRRHIcFa2IZFiyLYm33m3Ae25bqROThyUDt1KwZEWrY4nkSCpaEcmQiEun6DE2lFUBF3kmpjwTx+zQTE8it6EpGEUk3Q7uWkWt9x9gnd9Fpvh2Zcr7B1WyInegEa2IpMuy2cPpdWAEvi4GP1WbQN0OA62OJOIQVLQiclt2ezIj32vJO0lrCY3xZkm/dZSo9IjVsUQchopWRG4p6up5eo8OZan/GXpHlGTqu7vwCshvdSwRh6JztCKSpqO711N7VBArfM/wsWsbvvnomEpW5C5oRCsiN1k1ZwTdfx+Om7vBmkrv07jza1ZHEnFYKloRSZVsS2LE6GaMNDdQOd6LpX1XU6pKfatjiTg0Fa2IAHD57J/0HPcwawIu0zuqNJ8P34ZPQAGrY4k4PJ2jFRG2/TSDkAnl+Z/vZab59eCbD4+qZEUyiUa0Ik7MNE0+m/QEr1yeRXHTlS11v6ZG875WxxLJVVS0Ik4qOuIiT4+pzVzf47SOyM+MQVvId385q2OJ5DoqWhEndGjHKjrObc8RvwTG2BsxZPxqXNzcrY4lkiupaEWczNwvB9Lv+Cf4uBmsrTiGxl1ftzqSSK6mohVxEvHREbw6sh6TffZTJ9qPec+up3iFWlbHEsn1VLQiTuDwrtV0nfU4+wLieSWxBu+99z/cvX2tjiXiFFS0IrmYaZp8O+15Xjg1BR93g+9Lv03rJ0ZYHUvEqahoRXKpqIgLPDemDrN8/6JhVB5mPrOG4hUfsjqWiNNR0YrkQrs2zaPbst4c801khK0eb3ywBldPL6tjiTglFa1ILmLa7XwyoSuvRS6ksOHChiofUa/jK1bHEnFqKlqRXOLyqSM8OaEB3wdcoF1kIaa/soH8JStYHUvE6WmuY5FcYN2iD6k6qSJrfC8w0as9S8efVcmK5BAa0Yo4sLjoCF4f24RPPH6jgs2D71t9Q/WG3a2OJSL/oKIVcVB7fllEz0W9OJgnnv/GVOL94RvwDshvdSwR+RcVrYiDSbYl8dH4TrwVvZwCri78+MA7NO853OpYInILKloRB3Li8DZ6T2vBxoAIOkQWZtrLP5E/qKLVsUTkNlS0Ig7ANE1mTXuRF05MxvSCrwP60Oft6Rguup5RJKdT0YrkcJdPHeGFCY8yP+AUdWL9+a7vMkpVa2R1LBFJJxWtSA62eMbrPHfgA6762Rnt2owhHyzH1cPT6lgikgEqWpEc6PK5v/jv+KbM9QujeqIXa9t8R5V6nayOJSJ3QUUrksMsnvkmz/3+Hle97YxIrs/QMT9oSTsRB6aiFckhLp/9k/9OaMZcv+NUT/BibfPpVGmsySdEHJ2KViQHWPLtUJ49+GHKKNZWn6FjVuLu42d1LBHJBCpaEQudD9vPSxNbMT/gFNUTvVnb8muqNOxqdSwRyUQqWhEL2O3JTJ/8NIPPfEOsr8lIowlDxi7H3cvH6mgikslUtCLZ7MjudfT/rjMbAyKoH5eHad3nUO6hVlbHEpEsoqIVySaJ8TG8/1F7RsWvxccTvvDvyVNvfYOLq/4aiuRm+hsukg1+Wf0lT695kYN5EugSVZxPnl9OkbIhVscSkWygohXJQlfPh/HmhDZM8T5ACRdXVpR6gzbDR1sdS0SykYpWJAvY7cl8O/U5XjvxJVe8TP6bUI1RQ77Hv0Bxq6OJSDZT0Ypksj2/LOKFhU/xS8A1Ho73ZXKLL6nWsJvVsUTEIipakUwSGX6GYRPa8anrb+RzN5ie5wn6vPkVLm7uVkcTEQupaEXukWm3M2vaiww6PpWL3naejSnPqAHLyHdfWaujiUgOoKIVuQe7/zePl5b2Z1PgNWom+PB9408Jbf6k1bFEJAdR0YrchfMn9vPm5M587X2Y/J4GUwOfoN9b03VPrIjcRL8VRDIgITaKjyf1YHTk98R7wisJIbw9cDEBRUpaHU1EcigVrUg6mHY7S2a+yeB9H3HMP4m2MQX5qOe3PFirpdXRRCSHU9GK3MHuDXN4ZdnzbAiMIDjZgzVlR/FotzfAMKyOJiIOQEUrcgvHD27hrelPMNv/OPk9DT7z7Uz/od/i5ultdTQRcSAqWpF/CT93jNGTu/IZO3HxgqFJDzHk5XkEFtZ5WBHJOBWtyHVxUVf55NNevBe1iih3k75RZXi330xKVKxtdTQRcWAqWnF6tsR4vp3yLMNPz+SMbzJtYgoytuNkKtXvZHU0EckFVLTitJJtScz95lXePTSFo3mSqJXky6zKI2nQ4WWro4lILqKiFadjtyezZOZbDNs7gYN5Eqhs92Rp0FDaPTEaw8XF6ngiksuoaMVpmHY7K+eOYNjOD9kdEEs50525hV+k85vjNfG/iGQZFa3keqbdzqr5oxmx7QO2BUZT2nDj27xP0mPoZ7pVR0SynIpWci27PZlls4cx6reP+S0glvvdXJkW0Iu+r03B3dvX6ngi4iRUtJLrJNuSWPDta4w+OIX9eeIpY7jxVWBver32GR7eflbHExEnc8eiNQxjOtAGuGiaZqWsjyRyd5LiY5n91UuMOTaDP/IkUsH0YGb+p+k69BMdIhYRy6RnRPsN8CkwI2ujiNyda1fO8cWXz/Hxpe857ZdMVbsXC4oNpMObH+giJxGx3B2L1jTNjYZhBGV9FJGMOXtsLxO/fZ4pSVuJ9DRpmBTA1AdfpmW3t3WbjojkGDpHKw7n4K/L+WjxYL7z+oNkAzrGFmNwvRHUbPkfq6OJiNwk04rWMIz+QH+A+++/P7N2KwKkXEG8av5oJm6byJrAcLw9oH98MC93mUCZ0EetjicickuZVrSmaU4DpgGEhoaambVfcW7Xws/yzbcDmXR2KX/6J1HM3YVRRhOeeeYzCtxXzup4IiJ3pEPHkiP9sWsNny0awtfsIcoTaif5MrLQ03Ts+4HugRURh5Ke23vmAA2BAoZhnAaGm6b5VVYHE+eTEB/DkjnDmLbva34OvIq7G3SJKcmA+m9Tq/lTYBhWRxQRybD0XHXcPTuCiPP64/cNfLHwDb5J2MZlbztBrq6MNpryZJ8PKVqmmtXxRETuiQ4diyUSYqNYMnc40/Z/y88BV3AFHosrTP/gZ3m02xu4uHtYHVFEJFOoaCXbmKbJ7k3zmfHj+8yy7+Gyt0mQiyujkxvyZLf3KVqxltURRUQynYpWstyZ4/uYteBtZlxYw4E88Xi4Qtv44jxdpT+PdnldszeJSK6mopUsER1xkSXz32XGkfms97+MacAjNl8+9+5Mlx6jyVfiQasjiohkCxWtZJr42GusWTKOBb99xxKvMGI8oBSuvJ30ML3avMmDD7e2OqKISLZT0co9iY+LYs2y8Sz47TuWu/3FNU/I62HQI7EcvUOfo07bFzDc9GMmIs5LvwElwxKiI1m9dBwL9s5hufv1cnWBTgll6BL8BI07DNKkEiIi16loJV0unTnKyhXjWHF0JWu8zhDtAXndjJRyrdCLxh1exd3X3+qYIiI5jopW0mSaJgd3/sCK9Z+z/OImfs1zDdOAYi4u9EyqwONVetKk3UCNXEVE7kBFK6muXTnHhlWfs2b/Un5IOsRxfxsAIXYvhiXXoV29p6nepBeGq6vFSUVEHIeK1onZbIns/N8c1v46izXh2/nVLxKbK3i7QqP4grzm25g2LV+iRPDDVkcVEXFYKlonYrcn8/uvy9i4dS4bzmzhJ89zRHiZGCaEJHsxyFaTZpU78UiL/nj6B1odV0QkV1DR5iCmaRKVYONqTCJXYhK5GpvIlZgkrsYkEpVgIzbBRmxSMnGJycQk2IhLSiY2MZkEWzJ2O9hN8/pHyudu9gTyhm/EjNjIWeMQ+wOvEOmVslTwfYYrjSNLUtmjHuVq9eG+slUpGuhF0QAvPD30YyEikln0GzWbxCUmc/5aPOcj4zl/LY5zkfFciIxP+fNayp9XYhKx2c1b7sPHwxUfD1e8PVzx9XDD+/rXAd7uuBgGLlHH4PQaIq/t5CTH+T3vNWK8AW8oc9WVOheLEehWjSsBTTjpU55d7nZ2JQCb4mHTttTXyePlRtEAb4oFelGmoB9lCvnxQCE/HijoR15fTfYvIpIRKtpMFBWfRNjlWI6HxxB2OeXjeHgMJ8JjuRKTeNP2fxda4QAvyhXxp4CfJ/l8Pcjr45Hyp68H+Xw8CPR1x8/DDReX/1+PNTL8DPt+WcqOg2v59dJuthlnOelnA29w84Cq17zpm1yZBkUbUa9Bb4qUq3HT6ycl24mKtxERm8ilqATOXy/8cxEp/xA4dTWOrcfCiU+yp/43+Xw9CC6Wh2r3BaZ+5PfzzJpvqIhILqCizaCYBBvHL6eUZ1h4DMevF2pYeAyXo28s06IBXgTl96V5cGFK5PWhSJ6UQ7NFrn/4pOMQrWm3c+LPnfy8cyV7jv/C3quH2WNc4LhfUuo2JXGltq0IL7lW5aHKLQip3w3vvAXvuG93Vxfy+aaUeumCfmluY7ebnImI489L0fx1MZqjF6LZdyaSz37+k78H3yXyehNaMi91HyxIvQcLUDiP1x1fW0TEWRimeetDlXcrNDTU3LlzZ6bvN7vEJSYTFv7/I9KU0WnKSPVSVMIN2xbO40nJ/L6Uyu9LUAFfShXwIaiALyXz+eLtkf7bYOz2ZE4f28PhfT9x5PgODl86xP6EU+z1vkakZ8r/I8OEByPdqGYWomreClR7oC4hD3egSOkqmfr+0yM20cb+M9fYc+oqe05FsP34ldR/aJQr7E+9BwtQv2xBapfOj4ebS7bnExHJToZh7DJNMzTN55y1aOOTkjkRHpsyIg2P4UTq6DSW89fib9i2gJ9nSoGmlqnv9c990jUq/Ztpt3Ph9GHCDv/K8bDdHDm3nyNRYRw2L3LEO5a4f6wWFxAPFeP8qOpegmqFqlK1XH0qP9QW38L3Zda3IFPZ7SaHz0ex6eglNh29zPawKyTa7OTxcuPRikVoVbkIdR8sgKeb7sEVkdzHaYs2NtHGySuxnAiP5WT4jedOz12L559vPb+vB0HXC/TvUWlQfl9K5vfB3yt966XGx0Vx7vg+zobt58zZw4RdOkrYtZOEJVzguEskYd4JxP+jlw0TSkW5Us4WSHnP4pQvWIFypWpSvkojCj1QDcPFcUeCcYnJbPnzMqv2n2ftwfNci7cR6OPOY1WL0Tn0PoKL5cEwjDvvSETEAThu0c6aBW++CSdPwv33w+jR0LNn6tN2u0l4TCKnrqYU6YnwWE5ciUn5/ErsTYd58/q4pxzm/ceItFQBX0rm9yXA++YyTUpK4Mqlk1y5cJzwiycIDz9N+NUzXIm6wOXoS5xLuMzZ5AjOGjGc80jgqtfN38u8cVAq3osgAinlWYSgvKUoVSyYoAdqUKZyA7z889779ymHS7TZ2fLnZRbvPsPqA+dJtNmpUDQPfR4uyWPVimfoELuISE7kmEU7axb07w+xsakPJbq7M+fxtqwrW5GIOBuRcUkk2+wYGBgkY5jJBHiZBHpDoJeJvyf4udvxcbPh5ZqMYU8kLjGG6KTrH7YYom1xRCfHEW2PJ5oEoknimksS4R7JXPO89ffGPRmKxrpSNMmTYoY/Rd3zUcynCEUDS1CsyAMUK1GRkg/WJKDw/aCRW6rI2CSW7zvL7G0nOXTuGoE+7nSteR9PPlKKIgG6iEpEHJNjFm1QEJw4cdPDYQFQ6uV72zWAhw38ksDP5oJfsit+dnf88MDP8MTf1Zv8bv7k8wwkv08B8vsXIn9gMfLnK06+wkHkL/YA/oXuc+hDu1YzTZPtx6/wzS9hrD5wHjcXFzqHluDZBmW4L5+P1fFERDLEMYvWxQXSyGYCP377FiaAaWIaRsqI0cXA3dUDdzcPPNy8cHf3xMPdC3dPHzy8fHH39Mbdxw8f37z4BhTAwzePRpo5xMnwWKZs/IsFO09hmtA5tAQDm5bVbUIi4jAcs2hvMaKlZEkIC7u3fUuOdC4yjqn/O8asbSdwdTHoV7c0zzQone6L0URErHK7os25xz5Hjwaffx1C9PFJeVxypaIB3rzTLpj1rzTk0YpF+PTnP2k07n8s3X2GrPgHoYhIdsi5RduzJ0ybljKCNYyUP6dNu+GqY8md7s/vw6Tu1Vn2Qh2KB3oxcN4eun/xK39ejLI6mohIhuXcQ8ciQLLdZO6Ok3zw4xHiEpN5pVlZnq5XGlcXnV8XkZzDMQ8diwCuLgY9HyrJ+lcb0Lh8Id5bdZguU7dy/HKM1dFERNJFRSsOoYCfJ5/3CuHjrtU4eiGK1hM3sWzPGatjiYjckYpWHIZhGDxevTirX65PcLE8vDR3D28s+Z34pGSro4mI3JKKVhxO0QBv5jxdm2cblGH2tpN0nrKVC/9aCEJEJKdQ0YpDcnN1YWjL8nzRO5Rjl6Jp9+lmfj8daXUsEZGbqGjFoT1asTALn3skZQrHqb+w7uAFqyOJiNxARSsOr0LRPCx9oQ7lCvvzzMxdLP7ttNWRRERSqWglVyjo78msp2tTu3Q+Xpm/l6+3HLc6kogIoKKVXMTP043pfWvSIrgI7644yPTNKlsRsZ6KVnIVTzdXJvWoTovgIoz4/iDfbQ2zOpKIODkVreQ67q4uTOxenaYVCvH2sgMs3KVztiJiHRWt5Eoebi581jOEOg/kZ+iifWz845LVkUTESaloJdfydHPl8141eKCQH8/P+o0DZ3WfrYhkPxWt5Gp5vNz5+sma+Hu50e/bnVyOTrA6kog4GRWt5HpFA7z5oncoV2ISeWHWbyQl262OJCJOREUrTqFS8QDGdqjMtuNXGPvDYavjiIgTUdGK0+gQUoK+jwQxfctxTdUoItlGRStO5fVW5alYNA+vLdrHRa34IyLZQEUrTsXTzZWJ3asRm2jj1QV7MU3T6kgiksupaMXpPFDInzdbV2TT0cvM33nK6jgiksupaMUp9ax1P7VK5WPUykM6hCwiWUpFK07JxcXgvQ6VSbDZGb78gNVxRCQXU9GK0ypd0I8BjR9g1f7zbPnzstVxRCSXUtGKU+tXrzT35fNmxIqD2DSRhYhkARWtODUvd1febFWBIxeimLNDF0aJSOZT0YrTax5chIdK5eOTdUeJS0y2Oo6I5DIqWnF6hmEwqHk5Lkcn8N2vYVbHEZFcRkUrAtQMykf9sgX5fMNfRCfYrI4jIrmIilbkulceLcvV2CTmbj9pdRQRyUVUtCLXVbsvkIdK5WP65uNaSk9EMo2KVuQf+tcvzdnIeH74/ZzVUUQkl1DRivxDo3KFKF3Ql+lbwqyOIiK5hIpW5B9cXAx6PVSSvaciOHj2mtVxRCQXUNGK/EuHkOJ4uLkwd4cuihKRe6eiFfmXQB8PWlUqwpLfzhCfpAksROTeqGhF0tCxRgmiEmxsOHLJ6igi4uBUtCJpeLh0fvL7erBi31mro4iIg1PRiqTBzdWFVpWLsv7QBWI0U5SI3AMVrcgttKpclPgkO5uOaq1aEbl7KlqRWwgNyou/pxs/H75odRQRcWAqWpFbcHd1oV7ZAvx85CKmaVodR0QclIpW5DYalSvExagEDp+PsjqKiDgoFa3IbTxcJj8A249fsTiJiDgqFa3IbZTI60PxQG+2HQ+3OoqIOCgVrcgd1CqVj+3Hr+g8rYjcFRWtyB2E3B/I5ehEzkbGWx1FRByQilbkDoKLBwBw4EykxUlExBGpaEXuoEKRPLgYsF/L5onIXVDRityBt4crpQv6ceicilZEMk5FK5IOZQr6EnY5xuoYIuKAVLQi6RBUwJcT4bEk23XlsYhkjIpWJB1KF/AlMdnOmatxVkcREQejohVJhxJ5fQA4G6miFZGMUdGKpEPhPJ4AXLime2lFJGNUtCLpUCiPFwAXryVYnEREHI2KViQd/D3d8HRz4WKURrQikjEqWpF0MAyDPN7uRCfYrI4iIg5GRSuSTv6eblyLV9GKSMakq2gNw2hhGMYRwzD+NAxjaFaHEsmJ/L3ciFLRikgG3bFoDcNwBT4DWgIVge6GYVTM6mAiOY2nmyuJtmSrY4iIg0nPiLYW8KdpmsdM00wE5gKPZW0skZzH1cXAbrc6hYg4mvQUbXHg1D++Pn39MRGn4upiYFPTikgGpadojTQeu2nCV8Mw+huGsdMwjJ2XLl2692QiOYxhgKY6FpGMSk/Rngbu+8fXJYCz/97INM1ppmmGmqYZWrBgwczKJ5JjJNjseLrpQn0RyZj0/NbYATxoGEYpwzA8gG7A8qyNJZLzJNjseLm7Wh1DRByM2502ME3TZhjGi8BqwBWYbprmgSxPJpLDJCQl4+HvaXUMEXEwdyxaANM0fwB+yOIsIjlaVLwNf690/ZUREUmlE04i6XQlJpF8Ph5WxxARB6OiFUmH+KRk4pKSyeurohWRjFHRiqRDeEwiAPlUtCKSQSpakXQ4czUOgOKB3hYnERFHo6IVSYczEbEAFM+rohWRjFHRiqTD6Ssa0YrI3VHRiqTDX5eiKRrgpQkrRCTDVLQi6fDHhWgeLOxvdQwRcUAqWpE7SLab/HUpmrKF/KyOIiIOSEUrcgfHL0eTYLNTtohGtCKScSpakTvYfTICgGr3BVqcREQckYpW5A72nIrA39ONBwrq0LGIZJyKVuQOdp+MoMp9Abi4GFZHEREHpKIVuY2rMYkcOn+NWkH5rY4iIg5KRStyG1uPhWOaUPdBFa2I3B0VrchtbPnzMr4erlQpoQuhROTuqGhFbsE0TTYcucTDZQrg7qq/KiJyd/TbQ+QW9p+5xpmIOJoFF7Y6iog4MBWtyC2sOXgeFwOaVlDRisjdU9GKpME0TVb+fo6aQfm02LuI3BMVrUga9p2O5NilGB6vXtzqKCLi4FS0ImlY/NtpPNxcaF2lqNVRRMTB/V979xdjxV2Gcfz7dpcKRRAMiBWwUFtRJGK1NQiGNNYqKgETY9RYxeqVf1tj1FYTU2+M0abRRK0xVUuU1BikaTGKINp4oW1URFoKLW2VAi5dtKGQqkXg9eIMyboLVz2/+cUz30+y2TOzyZnn3XP2PDszZ2ctWmmcZ06eYvOuEd60eA7TIe49LwAABYlJREFUJ0+qHUfS/zmLVhrnZ/eP8OTTJ3jXFfNrR5E0ACxaaZz1v93PxbOn8vpLZtWOImkAWLTSGH8+cJSdB47y/mUXEeE/EZD07Fm00hi33vMo0yYP847XzKsdRdKAsGilxr4njrNl92E+sHwB03wTlKQ+sWilxrfueZQpk4a4dsXC2lEkDRCLVgIeOnycu3Ye4pplL/ZKUJL6yqKVgK9s2cvU5wzzkSsvqR1F0oCxaNV59z32D7bvHeXDV76Eme7NSuozi1addvLUaW7a/CAXPm8y1y733Kyk/rNo1Wk/uHc/e0aO8YXVi5ly/lDtOJIGkEWrzho99m9u2fowK186m1VLXlg7jqQBZdGqkzKTz915PydOneaLa17hVaAkFWPRqpM27TjEL/eM8uk3L2LhrKm140gaYBatOmfkqX9x0+bdXLFgphenkFScRavOGT32DHOmT+bmdy5l6DwPGUsqa7h2AKltS+fPYOv1KznPkpXUAvdo1UmWrKS2WLSSJBVk0UqSVJBFK0lSQRatJEkFWbSSJBVk0UqSVJBFK0lSQRatJEkFWbSSJBVk0UqSVJBFK0lSQRatJEkFWbSSJBVk0UqSVJBFK0lSQRatJEkFWbSSJBUUmdn/O404Auzv+x23axbw99ohKuny7NDt+Z29m5z92bsoM2ef7QtFinYQRMQfMvPy2jlq6PLs0O35nd3Zu6aN2T10LElSQRatJEkFWbTn9p3aASrq8uzQ7fmdvZucvSDP0UqSVJB7tJIkFWTRnkVErIqIhyLikYi4oXaetkTE/Ij4dUTsiYjdEXFd7Uxti4ihiPhTRPy0dpY2RcSMiNgYEXubx/91tTO1JSI+2TzfH4iIOyJicu1MJUXE9yJiNCIeGLPu+RGxLSL2NZ9n1sxYyjlm/2rzvN8VEXdGxIx+b9eiHScihoBvAm8BFgPviYjFdVO15iTwqcx8ObAM+GiHZj/jOmBP7RAVfB3YkpkvA5bSke9BRMwFPgFcnplLgCHg3XVTFXc7sGrcuhuA7Zl5KbC9WR5EtzNx9m3Aksx8JfAwcGO/N2rRTvRa4JHMfCwzTwA/AtZWztSKzBzJzB3N7eP0Xmzn1k3VnoiYB7wNuK12ljZFxHRgJfBdgMw8kZlH66Zq1TAwJSKGgQuAv1XOU1Rm/gZ4ctzqtcD65vZ64O2thmrJ2WbPzK2ZebJZvBeY1+/tWrQTzQUOjFk+SIfK5oyIWABcBtxXN0mrvgZ8BjhdO0jLLgaOAN9vDpvfFhFTa4dqQ2YeAm4GHgdGgKcyc2vdVFXMycwR6P3CDbygcp5aPgj8vN93atFOFGdZ16m3ZkfEc4GfANdn5rHaedoQEauB0cz8Y+0sFQwDrwZuzczLgKcZ3EOH/6M5F7kWWAi8CJgaEdfUTaUaIuLz9E6fbej3fVu0Ex0E5o9ZnseAH0oaKyIm0SvZDZm5qXaeFq0A1kTEX+mdLnhDRPywbqTWHAQOZuaZoxcb6RVvF7wR+EtmHsnM/wCbgOWVM9XwRERcCNB8Hq2cp1URsQ5YDbw3C/zNq0U70e+BSyNiYUScT++NEXdXztSKiAh65+n2ZOYttfO0KTNvzMx5mbmA3mP+q8zsxJ5NZh4GDkTEombVVcCDFSO16XFgWURc0Dz/r6IjbwQb525gXXN7HXBXxSytiohVwGeBNZn5zxLbsGjHaU6Kfwz4Bb0fuB9n5u66qVqzAngfvb25nc3HW2uHUis+DmyIiF3Aq4AvVc7TimYvfiOwA7if3mviQF8lKSLuAH4HLIqIgxHxIeDLwNURsQ+4ulkeOOeY/RvANGBb85r37b5v1ytDSZJUjnu0kiQVZNFKklSQRStJUkEWrSRJBVm0kiQVZNFKklSQRStJUkEWrSRJBf0XKm/QFZJtseIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "x2, v2 = pyaudi_integration2.y.T\n",
    "plt.plot([it.constant_cf for it in x],[it.constant_cf for it in v],label='original')\n",
    "plt.plot([it.constant_cf for it in x2],[it.constant_cf for it in v2],'r',label='simulation')\n",
    "plt.plot([it.evaluate(dweights) for it in x],[it.evaluate(dweights) for it in v],'g',label='differential learning')\n",
    "plt.plot(x[0].constant_cf, v[0].constant_cf, 'ro')\n",
    "plt.legend(loc=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Differential learning xf: \t11.843378895251684\n",
      "Real xf: \t\t\t11.924939296996673\n",
      "Mean Absolute Difference xf:\t0.08156040174498891\n",
      "\n",
      "Differential learning vf: \t3.5577687770021136\n",
      "Real vf: \t\t\t3.7620227657073806\n",
      "Mean Absolute Difference vf:\t0.204253988705267\n"
     ]
    }
   ],
   "source": [
    "print(\"Differential learning xf: \\t{}\".format(x[-1].evaluate(dweights)))\n",
    "print(\"Real xf: \\t\\t\\t{}\".format(x2[-1].constant_cf))\n",
    "print(\"Mean Absolute Difference xf:\\t{}\".format(D.max(D.abs(D.to_float(x2[-1]) - D.to_float(x[-1].evaluate(dweights))))))\n",
    "print()\n",
    "print(\"Differential learning vf: \\t{}\".format(v[-1].evaluate(dweights)))\n",
    "print(\"Real vf: \\t\\t\\t{}\".format(v2[-1].constant_cf))\n",
    "print(\"Mean Absolute Difference vf:\\t{}\".format(D.max(D.abs(D.to_float(v2[-1]) - D.to_float(v[-1].evaluate(dweights))))))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:python3.7]",
   "language": "python",
   "name": "conda-env-python3.7-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}